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Abstract. We study numerically the statistics of Poincare recurrences for the Chirikov standard map and 
the separatrix map at parameters with a critical golden invariant curve. The properties of recurrences are 
analyzed with the help of a generalized Ulam method. This method allows to construct the corresponding 
Ulam matrix whose spectrum and eigenstates are analyzed by the powerful Arnoldi method. We also 
develop a new survival Monte Carlo method which allows us to study recurrences on times changing by ten 
orders of magnitude. We show that the recurrences at long times are determined by trajectory sticking in a 
vicinity of the critical golden curve and secondary resonance structures. The values of Poincare exponents 
of recurrences are determined for the two maps studied. We also discuss the localization properties of 
eigenstates of the Ulam matrix and their relation with the Poincare recurrences. 

PACS. 05. 45. Ac Low-dimensional chaos - 05.45.Pq Numerical simulations of chaotic systems - 05.45.Fb 
Random walks and Levy flights 



1 Introduction 

The interest to understanding of transition from dynam- 
ical to statistical description of motion had started from 
the dispute between Loschmidt and Boltzmann, which 
is now known as the Loschmidt paradox pj[2]. The two- 
dimensional (2D) symplectic maps represent an excellent 
laboratory for investigation of how statistical laws appear 
in dynamical, fully deterministic systems. Their proper- 
ties have been studied in great detail during last decades 
both on mathematical (see e.g. [3l|4] and Refs. therein) 
and physical (see e.g. [5]|6l|7] and Refs. therein) levels of 
rigor. The case of completely chaotic behavior, appearing 
e.g. in Anosov systems, is now well understood |3114] but a 
generic case of maps with divided phase space, where is- 
lands of stability are surrounded by chaotic components, 
still preserves its puzzles. A typical example of such a 
map is the Chirikov standard map [5 ,6 which often gives 
a local description of dynamical chaos in other dynamical 
maps and describes a variety of physical systems (see e.g. 
[8]). This map has the form: 

y = y -\ sin(27rx) , x = x -\- y (modi). (1) 

27r 

Here x, y are canonical conjugated variables of generalized 
phase and action, bars mark the variables after one map 
iteration and we consider the dynamics to be periodic on 
a torus so that 0<x<1^0<y<l. The dynamics is 
characterized by one dimensionless chaos parameter K. 



For small values of K the phase space is covered by in- 
variant Kolmogorov-Arnold-Moser (KAM) curves which 
restrict dynamics in the action variable y. For K > Kg 
the last invariant golden curve with the rotation number 
r = Tg = {{xt - xo)/t) = {Vb- l)/2 is destroyed 0|1O] 
and it is believed that for K > Kg the dynamics in y 
becomes unbounded [TTJ[T2]. A renormalization technique 
developed by Greene and MacKay [QtlTO] allowed to de- 
termine Kg = 0.971635406 with enormous precision. The 
properties of the critical golden curve on small scales are 
universal for all critical curves with the golden tail of the 
continuous fraction expansion of r for all smooth 2D sym- 
plectic maps [To]. Here and below the time t is measured 
in number of map iterations (due to symmetry there is 
also a symmetric critical curve at r = 1 — at i^T^). For 
K > Kg the golden KAM curve is replaced by a cantori 
[13 which can significantly affect the diffusive transport 
through the chaotic part of the phase space p!4[ [T5]. At 
any K there are some chaotic regions in the phase space 
bounded by internal or isolating (at K < Kg) invariant 
curves. 

The dynamics inside a chaotic component of the phase 
space (x, y) is characterized by correlation functions whose 
decay ensures a transition from dynamical to statistical 
description. The decay of correlations is related to the 
probability to stay in a given region of phase space since 
for a trajectory remaining in a small region the dynam- 
ical variables are strongly correlated. This probability in 
its own turn is related to the statistics of Poincare re- 
currences. Indeed, according to the Poincare recurrence 
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theorem [16] a dynamical trajectory of a Hamiltonian sys- 
tem with bounded phase space always returns, after a cer- 
tain time, to a close vicinity of an initial state. However, 
the statistics of these recurrences depends on dynamical 
properties of the system. For a fully chaotic phase space 
a probability to stay in a certain part of a phase space 
decays exponentially with time being similar to a random 
coin flipping [3^4, • However, in dynamical maps with di- 
vided phase space, like the Chirikov standard map, the 
decay of probability of Poincare recurrences P{t) is char- 
acterized by a power law decay P{t) oc has /3 ^ 1.5 
whose properties still remain poorly understood. 

One of the first studies of Poincare recurrences in dy- 
namical Hamiltonian systems with two degrees of freedom 
was done in [T^ where an algebraic decay with an ex- 
ponent /3 = 1/2 was found. This exponent corresponds 
to an unlimited diffusion on an infinite one-dimesional 
line which is in contrast to a bounded phase space. This 
strange observation was explained in [18,19 as a diffu- 
sion in a chaotic separatrix layer of a nonlinear resonance 
which takes place on relatively short diffusion times. On 
larger times, which were not accessible to the computa- 
tions presented in [17 , this diffusion becomes bounded by 
a finite width of the separatrix layer and a universal alge- 
braic decay takes place with the exponent f3 ^ 1.5 corre- 
sponding to a finite chaos measure [18, 19[ . This algebraic 
decay of P{t) has been confirmed by various groups in 
various Hamiltonian systems [2Qll2T] . [22] . [24|l25] . [2S1I27] . 
[28l[29],[30l[3T]. 

One can argue that such a slow algebraic decay with 
/3 ~ 1.5 appears due to trajectory sticking near stable is- 
lands and critical invariant curves and leads to an even 
slower correlation function decay C{t) ~ tP{t) with a di- 
vergence of certain second moments. A sticking in a vicin- 
ity of the critical golden curve [10 is expected to give 
/3 ~ 3 [ 24( 1251 . being significantly larger than the average 
value 13 ~ 1.5. A certain numerical evidence is presented in 
[27j showing that long time sticking orbits can be trapped 
not only in a vicinity of a critical golden curve but also in 
internal chaotic layers of secondary resonances. 

Theoretical attempts to describe trapping in secondary 
resonances as renormalization dynamics on some Cayley 
type tree was started in [22 with recent extensions done 
in [28] [32 ] l33] . However, a detailed understanding of the 
intriguing features of Poincare recurrences in the Chirikov 
standard map and other similar maps is still missing. 

In this work we use a generalized Ulam method devel- 
oped in [34","35^ and combine it with a new survival Monte 
Carlo method trying to reach larger time scales and to 
obtain a better understanding of statistics of Poincare re- 
currences in the Chirikov standard map and the separatrix 
map. 

The paper is composed as follows: in Section 2 we con- 
struct the Ulam matrix based on the generalized Ulam 
method and study the properties of its spectrum, eigen- 
states and corresponding time evolution for the case of 
the Chirikov standard map. The survival Monte Carlo 
method is introduced in Section 3 and the properties of 
the Poincare recurrences are studied with its help compar- 



ing results with the Ulam method. In Section 4 we apply 
the above methods to the separatrix map and in Section 5 
the localization properties of the eigenstates of the Ulam 
matrix are analyzed. The discussion of the results is pre- 
sented in Section 6. 

2 Generalized Ulam method with absorption 

The Ulam method was proposed in 1960 [36|. In the orig- 
inal version of this method a 2D phase space is divided in 
Nd = M X M cells and ric trajectories are propagated on 
one map iteration from each cell j. Then the matrix Sij 
is defined by the relation Sij = Uij/ric where Uij is the 
number of trajectories arriving from a cell j to a cell i. By 
construction we have Sij — 1 and hence the matrix 
Sij belongs to the class of the Perron-Frobenius opera- 
tors (see e.g. [37]). This Ulam matrix can be considered 
as a discrete Ulam approximate of the Perron-Frobenius 
operator (UPFO) of the continuous dynamics. 
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Fig. 1. (Color online) The left panel shows the eigenvalue 
spectrum \j for the projected case of the UPFO of map ([T]) at 
K ^ Kg in the complex plane for M = 280 and Nd = 16609 
by red/gray dots (projected matrix dimension Np — 15457). 
The green/gray curve represents the circle |A| = 1. The right 
panel shows the number Nc of eigenvalues, with modulus larger 
than Ac, versus Nd in a double logarithmic representation for 
Ac = 0.5 (crosses), Ac = 0.66 (stars), Ac = 0.8 (open squares) 
and Ac = 0.9 (open circles). The straight lines correspond to 
the power law fits Nc ^ Nd with exponents f = 0.971 di 0.006 
(Ac = 0.5), V = 0.919 =b 0.005 (Ac = 0.66), v = 0.832 ± 0.010 
(Ac = 0.8) and h = 0.821 ± 0.021 (Ac = 0.9). The fits are done 
for the data with Nc > 50, M > 35 and M < 400 (Ac = 0.5), 
M < 800 (Ac = 0.66), M < 1120 (Ac = 0.8), M < 1600 (Ac = 
0.9), since the Arnoldi method provides only a partial spectrum 
of the eigenvalues with largest modulus for large values of M. 

According to the Ulam conjecture [36^ the UPFO con- 
verges to the continuous limit at large M. Indeed, this 
conjecture was proven for ID homogeneously chaotic maps 
[38 . Various properties of the UPFO for ID and 2D maps 
are analyzed in [39^40 ,[4 1,42 J. Recent studies [44^45 demon- 
strated similarities between the UPFO, the corresponding 
to them Ulam networks and the properties of the Google 
matrix of the world wide web networks. It was shown that 
in maps with absorption or dissipation the spectrum of 
the UPFO is characterized by the fractal Weyl law [46 . 

The coarse-grained cell structure of the original Ulam 
method corresponds to an effective noise and in case of 
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Fig. 2. (Color online) Partial spectrum Xj for the projected 
case of the UPFO of the map Q at = i^^ for M = 1600. 
The left panel shows all eigenvalues obtained by the Arnoldi 
method with ua = 5000. The insert of the right panel shows 
the blue/black square of the first zoomed range of the left panel; 
the blue/black square here is the second zoomed range shown in 
the main figure of the right panel. The eigenvalue with largest 
modulus Ao = 0.99994672216 is indicated by an arrow. The 
green/gray curve represents in all cases the circle |A| = 1. 



a divided phase space the noise induces an artificial dif- 
fusion between chaotic and regular regions. In [34 this 
problem was solved by replacing the random initial points 
by a very long chaotic trajectory and the transitions be- 
tween cells are accumulated along the chaotic trajectory 
that keeps the invariant curves and stable islands even in 
presence of the effective noise. Furthermore, the matrix 
size is also reduced since only cells which are visited at 
least once by the trajectory are kept. Here we use this ap- 
proach for the analysis of the Poincare recurrences keeping 
the same notations as in [34 . In particular, as in [Sj, we 
exploit the parity symmetry x ^ 1 — x and y ^ 1 — y 
allowing to limit the effective phase space to < x < 1, 
< y < 0.5 and therefore reducing the number of cells 
at a given cell size by a factor of two. In x direction we 
use therefore M cells and in y direction M/2 cells with 
M e {25, 35, ... , 1120, 1600} and the intermediate val- 
ues are multiples of 25 or 35 by powers of 2. 

To study the Poincare recurrences withing the Ulam 
method we introduce absorption of all trajectories with 
y < Vcut =0.05. Thus we generate the matrix S using one 
trajectory iterated by the map up to the iteration time 
t = 10^^ (as in [34l; this corresponds to the closed system 
without absorption and we call this the symplectic case). 
After that the matrix size is simply reduced only to 
those cells with y > ycut that gives the projected matrix 
dimension Np and matrix Sp. The matrix size of this pro- 
jected case is smaller approximately by 7%. We find, for 
M < 1600, an approximate dependence Nd ^ 0.39M^/2 
and Np ^ 0.36M^/2. This corresponds to the usual esti- 
mate of the chaos measure being around 39% in agreement 
with the results of Chirikov [^ . For the symplectic case we 
have the maximal eigenvalue A = 1 while in the projected 
case with absorption we are getting |A| < 1. 

The spectrum Xj of the projected case with matrix 
Sp is shown in Fig. [l] The spectrum is obtained by the 
direct diagonalization of the matrix Sp that can be done 
numerically up to M = 280. It can be compared with the 



corresponding spectrum of the symplectic system shown in 
Fig. 2 of [34|. The global spectrum structure of S for the 
symplectic case is similar to the projected case. Indeed, 
the absorption is relatively weak and does not affect the 
global properties of motion. However, with absorption the 
measure is not conserved and the remaining non-escaping 
set forms a fractal set with the fractal dimension d < 2 
(see e.g. ffl|47!). 

In the case of Ulam networks on fractal chaotic re- 
pellers the spectrum of UPFO Sp is characterized by the 
fractal Weyl law with the number of states Nc in the ring 
Ac < |A| < 1 growing with the matrix size Nd as Nc oc N^ 
(here for simplicity we use the size A^^^ of the symplec- 
tic case, for the projected case we have simply to change 
Np ^ 0.93A^(i). It can be argued that the fractal dimen- 
sion do of the invariant repeller set determines the expo- 
nent V = do/ 2 [46]. Examples of dependencies Nc vs N^, 
are given in Fig. [l]for various values of Ac- Definitely we 
have u < 1 but there is an evident dependence on Ac with 
a decreasing value of at Ac ^ 1. We attribute this to 
the fact that at Ac ^ 1 we are dealing with long sticking 
trajectories whose measure decreases with time. 
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Fig. 3. (Color online) The left panel shows the density 
p{X) of eigenvalues, being normalized by f p{X)d^X = 1, 
of the UPFO for the map ^ Sit K = Kg in the com- 
plex plane as a function of the modulus |A| for M = 280 
for the symplectic case (upper curve, crosses) and the pro- 
jected case (lower green curve, stars). The other curves are 
partial (non- normalized) densities for the projected case and 
the values M = 400, 560, 800, 1120, 1600 and the number 
of used eigenvalues (obtained by the Arnoldi method) is 
UA = 12000, 8000, 8000, 6000, 5000 respectively. The right 
panel shows the decay rates 7^- = — 21n(|Aj|) versus level num- 
ber j for the UPFO eigenvalues Aj, with M — 1600 and 
Nd = 494964. The red/gray crosses correspond to the UPFO of 
symplectic case and the blue/black squares correspond to the 
projected case (data points for this case are shifted to one posi- 
tion to the right). The green curve corresponds to the quadratic 
dispersion law 7^ ~ 71 which is approximately valid for the 
diffusion modes with < j < 5 and where 71 is taken from the 
UPFO of the symplectic case. 

Here we should point out that the data for M > 400 
corresponding to A^^^ > 30000 are obtained from the Arnoldi 
method [48 which allows to find the eigenvalues for ma- 
trix sizes up to Nd ~ 10^. However, only a finite num- 
ber of eigenvalues with largest |A| can be determined nu- 
merically using Ua = 12000, 8000, 8000, 6000, 5000 (for 
M = 400, 560, 800, 1120, 1600 respectively and with ua 
being the used Arnoldi dimension). A more detailed de- 
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scription of the Arnoldi method for the UPFO is given 
in [34 . An example of the spectrum A obtained with the 
Arnoldi method at the largest value of M = 1600 is shown 
in Fig. [2] Here Nd = 494964 and Np = 458891. We 
find that the maximal eigenvalue for the projected case is 
Ao = 0.99994672216 corresponding to a slow escape rate 
at large times. As in [34 for the symplectic case without 
absorption, we obtain also for the case with absorption two 
type of eigenmodes: "diffusion modes" with real eigenval- 
ues close to 1 and whose eigenvectors are rather extended 
in phase space (with some decay for cells close to the ab- 
sorption border) and "resonant modes" with complex or 
real negative eigenvalues and which are quite well local- 
ized around a chain of stable islands close to an invariant 
curve. It turns out that many of the resonant modes (those 
"far" away from the absorption border), coincide numeri- 
cally very well with corresponding resonant modes for the 
case without absorption already found in [34 . 

The dependence of the density of eigenvalues p(|A|) on 
|A| is shown in Fig. [sj We see the proximity between the 
symplectic and projected cases not only in density p but 
also in a slow relaxation of the diffusion modes with relax- 
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Fig. 4. (Color online) The left panel shows the rescaled level 
number j/Nd versus the decay rate 7^, in a double logarithmic 
scale, for the map Q at Kg with M = 1600 and Nd = 494964. 
Red/lower data points correspond to the UPFO projected case 
and green/upper data points correspond to the UPFO sym- 
plectic case. The two straight lines correspond to the power 
law fits j/Nd ^ 0.052745 7^-^^°^ (symplectic case) and j/Nd ^ 
0.0415707^-^^^^ (projected case) for the data in the range 
0.04 < 7 < 0.3. The statistical error bound of the exponents 
obtained from the fits is close to 0.1% in both cases. The right 
panel shows the decay rates 'jj{M) for j = (red crosses), j = 5 
(green open squares) of the UPFO projected case in a double 
logarithmic scale. The lower/pink straight line corresponds to 
the power law fit 7o(M) ^ 0.72 M"^-^^ and the upper/light 



obtained for the range 400 < M < 1600). The black/curved 
line corresponds to the other fit 70 (M) — /(M) — i+b/m 
with D = 0.162, C = 165 and B = 17.0 (fit obtained for the 
range 25 < M < 1600). We mention that 75 corresponds for 
M > 400 to a resonant mode whose eigenvector is strongly 
localized close to the three stable islands of the resonance 1/3. 
However, for M < 280 75 corresponds to a different mode and 
the resonant mode at 1/3 is associated to 77 (M = 280), 713 
(M 200), 717 (M 140) and 723 (M 100) which are 
shown as four additional data points (blue stars). 



ation rates ^ 71 (7^- = -21n|A^ |) provided we iden- blue straight line to the fit 75 (M) ^ 389 M" • (both fits 
tify 7j+i of the symplectic case with of the projected 
case because 70 of the symplectic case is simply zero and 
the relaxation rate 71 to the ergodic state of the symplec- 
tic case corresponds roughly to the exponential long time 
escape rate 70 of the projected case. The proximity of the 
two cases is also well seen in the dependence of integrated 
density of states Pe{i) = j/^d on 7j shown in Fig. H (here 
j is a number of eigenvalues with 7 < 7^). In both cases we 
have the algebraic dependence Pe{i) ^ 7^ with f3 ~ 1.5. 
In [34] it was argued that this exponent is the same as 
for the exponent of decay of Poincare recurrences P{t). 
These data show that an introduction of small absorp- 
tion dit y < ycut does not produce significant modification 
for trajectories trapped for long times in a vicinity of the 
critical golden curve or other secondary islands located far 
away from the absorption band. 

The lowest eigenvalues such like 70 and 75 decrease 
algebraically with the increase of M as it is shown in right 
panel of Fig. [4j In the fit range 400 < M < 1600 we have 
a power law 70 (M) « 0.72 M~^-^^ but taking into account 
the curvature for the interval 25 < M < 1600 the modified 



fit 7o(M) = ^ \X%^iM with D = 0.162, C = 165 and B = 
17.0 seems to indicate a behavior 70 (M) cx in the 

limit M ^ 00. This behavior is similar to the one found 
in [34] for 71 in the symplectic case (where 70 is simply 0). 
On the other hand the resonant mode 75 obeys the power 
law 75(M) ^ 389 M-^-^^ which is valid for the interval 
100 < M < 1600 if we use for the smaller values of M 
not 75 but the resonant mode localized to the same chain 
of resonant islands which may have a different eigenvalue 
index (see Fig. |4] for details). The comparison of these 
decays indicate that eventually at very large values of M, 
far outside the range numerically accessible by the Arnoldi 
method, the resonant modes become dominant over the 
diffusion modes. The limit 7 ^ for M ^ 00 is related 
to long sticking trajectories near critical invariant curves 



which restrict the chaos component and whose phase space 
structure can be better resolved with decreasing cell size 
1/M. As in [34] we argue that these lowest modes are 
affected by the effective noise present in the Ulam method. 
Due to that we do not have a clear explanation for this 
algebraic decay. However, the fact that jj (at fixed value of 
j) vanishes with increasing M indicates that the limit texp 
in time, when the statistics of Poincare recurrences P{t) 
obtained from the UPFO becomes exponential, increases 
as well according to texp 0^ 7o^^ and therefore we expect 
to recover the power law decay of P{t) for M ^ 00 (see 
below in Section 3). 

With the help of the Arnoldi method we find certain 
eigenstates corresponding to eigenvalues of the matrix Sp 
and satisfying the equation 



(2) 



i=0 



Examples of two eigenmodes |?/^o| and |?/^29| are shown in 
Fig. [5] The state {tpol corresponds to the first diffusive 
mode mainly located in a vicinity of the critical golden 
curve while |?/^29| corresponds to the mode located near a 
resonant chain with rotation number r = 2/7. 
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Fig. 5. (Color online) Density plot of the modulus of the 
eigenvector components IV^I of the UPFO projected case of map 
([T]) at Kg with M — 1600 for the two modes with eigenval- 
ues Ao = 0.99994672 (left panel) and A29 = -0.22008951 + 
i 0.96448508 ^ IA29I e^^^^^^^^ (right panel). The density is 
shown by color with red/gray for maximum and blue/black 
for zero. 




Fig. 6. (Color online) Time dependent probability density 
calculated by ^(t) = (5'p)*V^(0)/ || (5'p)*V^(0) ||i where Sp 
is the UPFO for the projected case for M = 1600, '0(0) an 
initial vector with '0^(0) = 81^1^ and Iq being the index of the 
cell at xq — yo — 0.0625 and || • • • ||i is the 1-norm defined by 
II -0 II 1= The densities are shown for t = 40 (left panel) 

and t — 400 (right panel). In the limit t ^ 00 the vector '0(t) 
converges to the eigenvector of maximal eigenvalue Ao shown 
in the left panel of Fig. [5] The full convergence is achieved 
for t > 40000 so that for these times the density plot of ip(t) 
remains unchanged at the given color-resolution. 

It is also interesting to follow how the probability ini- 
tially placed in one cell io evolves with time. Of course, 
the total probability starts to decay due to absorption but 
by renormalizing the total probability back to unity after 
each map iteration we obtain its evolution in phase space. 
At large times we have convergence to the state i/jq with 
maximal Aq but at intermediate times we see the regions 
of phase space which contribute to long time sticking and 
long Poincare recurrences. Two snapshots are shown in 
Fig. [gI The videos of such an evolution for the maps ([T]) 
and O) are available at [35 . 



3 Poincare recurrences 

by survival Monte Carlo method 

The numerical computation of the Poincare recurrences 
counting the number of crossing of a given line (e.g. ^ = 0) 
in the phase space is known to be a very stable numerical 
method since the integrated probability of recurrences on 
a line at times larger than t is positively defined (see e.g. 
p!8 l fT9] . [24 | [28] ) . However, at large times the direct numer- 
ical computation becomes time consuming. 

With the aim to reach larger times we present here 
a new method to calculate the statistics of Poincare re- 
currences of map such as the Chirikov standard map (IT]). 



We will call this method the Survival Monte Carlo method 
(SMCM). The idea of this method is to chose a certain, 
quite large number A^^ 1, of initial conditions randomly 
chosen in some small cell close to an unstable fix point and 
to calculate in parallel the time evolution of these trajecto- 
ries. At the initial time t = we put the Poincare return 
probability to P(0) = 1 and the number of trajectories 
to N{0) = A^^. At each time tk, when a given trajectory 
escapes in the absorption region y < ycut = 0.05 of the 
phase space, we put P{tk + 1) = P{tk) {N{tk) - 1)/N{tk) 
and N{tk + 1) = N{tk) — 1, otherwise we simply keep 
P{tk + 1) = P{tk) and N{tk + 1) = N{tk). When the 
number of remaining trajectories N{tk) drops below a 
certain threshold value Nf (typically chosen such that 
A'i ^ A^j ^ 1) we reinject a new trajectory close to 
one of the other remaining trajectories with a small ran- 
dom deviation: Xnew(^) ^ [xi{t) — Xi{t) -\- e/2] and 
ynew{t) G [yi{t) — yi{t) -\- e/2]. The main idea is to 
keep a typical statistics of trajectories at a given time t 
and to concentrate the computational effort on the very 
long and rare trajectories without wasting resources on the 
more probable trajectories with short times of Poincare 
recurrences. 




Fig. 7. (Color online) Statistics of Poincare recurrences P(t) 
of the map ([T]) calculated by the SMCM as survival probability 
after times larger than t (data are shown in double logarith- 
mic scale). The number of initial trajectories is Ni — 10^ and 
the number of final trajectories is Nf — 100 (left panel) or 
Nf — 1000 (right panel). The initial positions are randomly 
chosen in a cell of size (1600)"^ X (1600)"^ at the position 
xo — yo — 0.0625, here the small random deviation for rein- 
jected trajectories is ^ £ = 10~^^. In both panels the results 
for P{t) are shown for 10 realizations with different random 
seeds. The horizontal dotted line indicates the limit proba- 
bility Nf/Ni = 10"^ (left panel) or Nf/Ni = 10"^ (right 
panel) below which the reinject ion of trajectories is applied. 
The two realizations in the left panel which drop below the 
shown range (of P{t) > 10~^^) "saturate" eventually at the 
values P(10^^) 2 x 10"^^ or P(10^^) ^ 10"^^ 



In this method the proper choice of e is important. 
On one hand e should not be too small in order to avoid 
too strong correlations between the trajectories and on 
the other hand it should be very small in order to avoid 
an uncontrolled too strong diffusion into regions too close 
to stable islands where the trajectories may be trapped 
stronger and longer as they should be without the random 
deviations. Fortunately in the chaotic region even a mod- 
est Lyapunov exponent ensures exponential separation of 
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trajectories and choosing a very small value of e one may 
hope to reduce the correlation between the injected trajec- 
tory and its reference trajectory after a modest number of 
iterations. Furthermore at longer times the average time 
between the escape of two trajectories becomes very large 
that helps to reduce these correlations. 
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Ulam, M = 


soo^-^^^s. 
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= 1600 -A. \ 




Nf= 1000 
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Fig. 8. (Color online) The left panel shows the average over 
10 random realizations of the statistics of Poincare recurrences 
P{t) of the map ([T]), obtained by the SMCM for survival prob- 
ability and shown in Fig.]?] The upper/red curve, for t < 10^°, 
corresponds to Nf — 1000. The next lower/green curve, for 
t < 10^^, corresponds to Nf — 100. The lowest blue curve, for 
t < 1.7 X 10^, corresponds to the data of Refs. |251l26j obtained 
by a direct computation of the statistics of Poincare recur- 
rences. The dashed straight line indicates a power law behavior 
P oc t~^'^ . The right panel compares the statistics of Poincare 
recurrences P(t) obtained by the SMCM for Nf = 1000 to 
P(t) obtained by the Ulam method for M = 400, 800, 1600. 
At large times t > texp 10^ — 10^ the curves obtained by the 
Ulam method show an exponential behavior P{t) ^ Xq deter- 
mined by the largest eigenvalue of the UPFO for the projected 
case. 

We have chosen the parameters e = 10~^^, Ni = 10^ 
and the two cases Nf = 100 and Nf = 1000. For Nf = 100 
we have been able to iterate up to times 10^^ and for 
Nf = 1000 up to times 10^°. We mention that at the 
"larger" value £ = 10~^^, we observe a significant effect 
of artificial saturation in P{t) (i. e. no escapes of trajec- 
tories) at longer times because the trajectories penetrate 
"too strongly" into the regions very close to the stable is- 
lands or the critical curve. When we choose e = 10~^^ this 
artificial saturation effect is strongly reduced. We calculate 
in parallel different realizations of P{t) with respect of the 
random variables (for the initial conditions, for the ran- 
dom deviations of the reinjected trajectories and for the 
random choice at which remaining trajectory the reinjec- 
tion happens). The comparison of obtained data shows 
that the distribution P{t) is stable at small and large 
times. But at very large times it turns out that the fluc- 
tuations become quite strong. 

Examples of the survival probability P{t) obtained for 
10 different realizations with Nf = 100 (left panel) and 
Nf = 1000 (right panel) are shown in Fig. ItJ Of course 
the fluctuations appear for Nf = 100 at snorter times 
{t - 10^ - 10^) as compared to Nf = 1000 {t - 10^ - 10^). 

We note that the SMCM allows us to determine the 
survival probability P{t). Its comparison with the statis- 
tics of Poincare recurrences computed by the usual method 



Fig. 9. (Color online) Density plots of the trajectories of 
the SMCM (with Nf = 1000) for the map ^ for various 
times t and random realizations. All density plots are obtained 
from a histogram of 10^ data points and using a resolution of 
800 X 400 cells for the phase space < x < 1 and < y < 0.5. 
The data points are obtained by iterating N{t) trajectories 
(with N(t) = P(t) Ni for P(t) > 10"^ and N{t) = Nf for 
P(t) < 10"^) from ttot + At with At = 10^ /N(t). The left 
four panels and the upper right panel correspond to one par- 
ticular random realization at t = 10^, 10^, 10^, 10^, 10^° and 
the three lower right panels correspond to three other ran- 
dom realizations at t = 10^*^. For short times t < 10^ there 
is no significant difference between the density plots for differ- 
ent random realizations at a given time. More detailed density 
plots for intermediate times and higher resolution figures are 
available at |35j. 

[T81I24]. [251 126] is shown in Fig.jsj We see that both meth- 
ods give the same behavior P{t) with a small shift in time 
related to different initial conditions. The equivalence of 
both methods is rather clear: in both methods the prob- 
ability is determined by long sticking trajectories; both 
methods consider the recurrences to the lines y = or 
y = 0.05 which are close to each other. 

The decay of P{t) averaged over 10 random realiza- 
tions is shown in Fig. [Sj In general we see that the SMCM 
allows to reach extremely long times with t = 10^^ for 
Nf = 100 and t = 10^° for Nf = 1000. For Nf = 100 
we see that the fluctuations start to be important for 
t > 10^ while the case with Nf = 1000 remains stable 
up to t = 10^^. This allows to obtain the behavior of P{t) 
for times being about one order of magnitude larger com- 
pared to previous numerical simulations. 

For the case Nf = 1000 in Fig. [8] the algebraic fit 
of data in the range 10^ < t < 10^° gives the Poincare 
exponent (3 = 1.587 ± 0.009. For Nf = 100 case we find 
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(3 = 1.710 ± 0.017 for the range 10^ < t < 10^^ The 
formal statistical error is rather small in both cases but 
it is clear that for Nf = 100 we start to have an effect 
of strong fluctuations due to long sticking around islands 
and thus the reliable value of P is given by the case with 
Nf = 1000. 

The survival probability P{t) can be also computed 
using the Ulam method at various sizes of discrete cells 
determined by M. The results obtained by the generalized 
Ulam method and by the SMCM are shown in the right 
panel of Fig. ^ The comparison shows that both methods 
give the same results but the SMCM is much more efficient 
allowing to follow the decay P{t) up to significantly larger 
times since for the Ulam method we expect the decay P{t) 
only to be accurate for t < texp 7o^^ because for t > texp 
it becomes exponential P{t) ex Aq = exp(— 7ot/2). The 
data of Fig. [S] clearly shows that texp increases with M in 
accordance with the decay of 70 obtained from Fig. |4] 

Using the SMCM we can follow the evolution of the 
survival probability as a function of time showing the den- 
sity plot of long sticking trajectories. Examples of such 
distributions are shown in Fig. [9j These Figs, show that 
at short times t < 100 the trajectories are not yet able to 
cross the cantori barriers and remain relatively far from 
the golden curve, at larger times t = 10^, 10^, 10^ the prob- 
ability becomes concentrated close to the golden curve. 
But at very larger times t = 10^^ we find trajectories 
sticking in a vicinity of the golden curve or other sec- 
ondary resonances. Thus we see that at long time P{t) 
has contributions not only from the vicinity of the critical 
golden curve but also from other secondary resonances. 
In this respect, our conclusion confirms a similar one ex- 
pressed in [27 obtained from simulations on shorter time 
scales. 



4 Separatrix map with critical golden curve 

To show that the previous case of the Chirikov standard 
map represents a generic situation we also study the UPFO 
of the projected case for the separatrix map [6 , defined 
by: 

^ = y + sin(27rx) , x = x -\ ln(|^|) (modi). (3) 

27r 

This map can be locally approximated by the Chirikov 
standard map by linearizing the logarithm near a certain 
yo that leads after rescaling to the map ([T]) with an effec- 
tive parameter i^eff = ^/|^o| [6 • As in [34J we study the 
map (|3| at ylc = 3.1819316 with the critical golden curve 
at the rotation number r = = (a/5 — l)/2 = 0.618... The 
construction of the matrix S is described in [34] , its size is 
given by an approximate relation ^ 0.78M^/2 for the 
phase space region 0<x<l,0<?/<4 (symplectic case 
and using the symmetry: x ^ x -\- 1/2 (mod 1), ^ — ?/). 
The absorption is done for y < ycut = 0.4 corresponding 
to 10% of the maximal possible value of y. Thus for the 
UPFO for the projected case we have Np ^ 0.68MV2. 
In fact we have 2{Nd — Np)/M'^ = 0.1 since all part of 



the phase space is chaotic at < < ycut and all cells in 
this region were occupied by the Ulam method. Thus for 
M = 1600 we have A^^ = 997045, Np = 869045. 




Fig. 10. (Color online) The left panel shows the rescaled level 
number j/Nd versus the decay rate 7^, in a double logarith- 
mic scale, for the separatrix map ^ at Ac with M = 1600 
and Nd = 997045. Red/lower data points correspond to the 
UPFO for the projected case and green/upper data points cor- 
respond to the symplectic case. For the symplectic case the 
data points are shifted up by a factor 2 to separate the two 
data sets. The two straight lines show the power law fits j/Nd ~ 
0.0141737^-^^^^ (symplectic case) and j/Nd ^ 0.0142077^-^°^^ 
(projected case) for the range 0.04 < 7 < 0.3. The statisti- 
cal error of the exponents is close to 0.2% in both cases. The 
right panel shows the decay of 7^- (M) with M for j = (red 
crosses), j — 2 (green open squares) for the UPFO for the 
projected case of map (|3|. The lower /blue straight line cor- 
responds to the power law fit 7o(M) ^ 2.26 M~^"^^ and the 
upper/pink straight line to the fit 72(M) ^ 1,95M"°-^^ (for 
the range 400 < M < 1600). The eigenvector corresponding to 
72 is localized near the two stable islands of the resonance 1/2. 



In Fig. [Toj in analogy to Fig. [4j we show the depen- 
dence of integrated number of eigenvalues j/Nd on 7^ = 
— 21n|Aj| for the symplectic and projected cases of the 
UPFO of the map (|3|. In both cases we have approxi- 
mately the same dependence with the algebraic exponent 
P ^ 1.5 which works for the range 0.04 < 7 < 0.3. The 
minimal values of 7 (e.g. 70 and 72) drop approximately 
inversely proportionally to M. As for symplectic case ^34] 
we attribute this decrease with M to a finite size coarse- 
graining effect of the Ulam method. As in [34 , we ar- 
gue that the exponent /3 for a more physical intermediate 
range of 7 is directly related to the Poincare exponent. 



Fig. 11. (Color online) Density plot of the modulus of 
the eigenvector components of the UPFO for the projected 
case of the map Q at M = 1600 for the two modes with 
Ao = 0.99972660 (left panel) and A77 = -0.49158775 + 
i 0.85153885 IA77I e^^^^^/^^ (right panel). 
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Examples of two eigenmodes at Aq and A77 are shown 
in Fig. 11 In the first case we have an eigenmode of diffu- 
sive type similar to Fig. |5] while in the latter case we have 
an eigenmode concentrated around unstable fix points of 
resonance 1/3 (see corresponding state of symplectic case 
in bottom left panel of Fig. 11 in [34]). 





10° 10^ 10^ 10^ 10^ 10^° 



Fig. 12. (Color online) The left panel shows the average over 
10 random realizations of the statistics of Poincare recurrences 
P(t) of the map ([3|, obtained by the SMCM. The red curve, 
for t < 10^°, corresponds to Nf = 1000. The green curve, for 
t < 10^^, corresponds to Nf = 100. The upper/blue curve, 
for t < 2.8 X 10^, corresponds to the data shown in 25 using 
a direct computation of the statistics of Poincare recurrences. 
The right panel compares P(t) SMCM data for Nf = 1000 (red 
curve in left and right panels) with P{t) obtained by the Ulam 
method for M = 400, 800, 1600. At large times t > texp ^ 
2x 10"^— 2x 10^ the Ulam method leads to an exponential 
decay P{t) ^ Xq determined by the largest eigenvalue of the 
UPFO for the projected case. 



The comparison of the statistics of Poincare recur- 
rences obtained from the map (p|by the SMCM and the 
usual method are shown in Fig. |12[ The data of the usual 
method obtained in [25] allows to follow the decay of P{t) 
up to t = 2 X 10^, while with the SMCM we reach times 
t = 10^0 with Nf = 1000 and t = 10^^ with Nf = 100. 
We have a good agreement between three curves for the 
range 100 < t < 10^ with a certain constant displacement 
in logio ^ of ^^^^ from the usual method compared to the 
SMCM data. This shift appears due to different initial 
conditions but apart of this shift all oscillations of P{t) 
curve are well reproduced. This shows that both methods 
works correctly. However, with the SMCM we are able to 
reach times being by one to two orders of magnitude larger 
than previously. 



The algebraic fit of SMCM data in Fig. [12] gives /3 = 
1.855 ± 0.004 for Nf = 100 (range 10"^ < t < 10^^) and 
p = 1.706 ± 0.004 for Nf = 1000 (range 10^ < t < 10^°). 
In both cases the statistical error is rather small but there 
are visible fluctuations which become to be significant at 
t > 10^ for Nf = 100 even if they are smaller compared to 
the similar case of map ([T]) shown in Fig. [sj Due to that 
one should take as the reliable value /3 = 1.706 that shows 
a noticeable difference from the value /3 = 1.587 found 
above for the Chirikov standard map at = Kg. 

The comparison of the SMCM data for P{t) with the 
results of the Ulam method are shown in the right panel 
of Fig. 12 As it was the case for the similar comparison 
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Fig. 13. (Color online) Density plots of the trajectories of 
the SMCM with Nf = 1000 for the map ^ for various times 
t and various realizations. All density plots are obtained by a 
histogram of 10^ data points with a resolution of 800 x 400 
cells for the phase space < x < 1 and < ^ < 4. The 
data points are obtained by iterating the N(t) trajectories 
(with N(t) = P(t) Ni for P(t) > 10"^ and N[t) = Nf for 
P(t) < 10"^) from ttot + At with At = 10^ /N(t). The left 
four panels and the upper right panel correspond to one par- 
ticular random realization at t = 10^, 10^, 10^, 10^, 10^° and 
the three lower right panels correspond to three other ran- 
dom realizations at t = 10^*^. For short times t < 10^ there is 
no significant difference between the density plots for different 
random realizations at a given time. 



results but the Ulam method works only for time scales 
being significantly smaller than those reached with the 
SMCM. 

Finally, as in Fig. [9] we show in Fig. [13] the density 
distribution obtained for various realizations and various 
times of the map ([3|. The situation is similar to Fig. [9j 
at short times the density is bounded by cantori barriers, 
at large times it reaches the critical golden curve and at 
even larger times we see that the density is located near 
the critical golden curve or other secondary resonances 
depending on the realization. 



shown in Fig. [S] we find that both methods give the same 



5 Properties of eigenstates of Ulam matrix 

Let us now try to analyze how the decay of Poincare recur- 
rences is related to the properties of the (right) eigenvec- 
tors il){x^y) of the UPFO for the projected case. For this 
we determine the x-average of the eigenvector amplitude 
around a given position xq over a band of 1% width of the 
whole x-range: (|V^(^)|) = 100 E|z^a.|<o.oo5 IV^(^o + 
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Fig. 14. (Color online) The localization properties in y- 
direction for certain eigenvectors of the UPFO for the projected 
case for the maps ([T]) (left column) and ^ (right column). The 
panels in the second and fourth row show the averaged mod- 
ulus {\ip{y)\) of the eigenvector components within a band of 
1% width of the whole x-range at a certain x = xq. The global 
structure of the corresponding eigenstates is shown in the cor- 
responding first and third panels (counting from the top; the 
red vertical thick line indicates the range of x-values where 
the average has been performed for each value, M = 1600). 
Data are shown for M = 400 (cyan/highest curve), M = 560 
(pink/second curve), M = 800 (blue/third curve), M = 1120 
(green/fourth curve) and M = 1600 (red/lowest curve). In the 
right panel of the second row the data for different values of 
M approximately coincide and only the data for M = 1600 are 
shown by a full (red) curve; other M values are shown as iso- 
lated data points for M = 1120 (green crosses), M = 800 (blue 
stars), M = 560 (pink squares) and M = 400 (cyan circles). 
For M = 1600 the eigenvectors, shown in the density plots of 
the first and third row, correspond to the modes A4 and A31 of 
the map ([T]) (left column) and to the modes A2 and A17 of the 
map (|3| (right column); for other M we show corresponding 
eigenvector located at the same resonances. 



Ax^y)\. The ?/-dependance of this average allows to vi- 
sualize the localization properties of the eigenstate in y- 
direction. In Fig 



14 we show this quantity for two ex- 
amples for each of the maps ([T]) and ^ and for different 
values of M between 400 and 1600. 

For the case of the map ([T]) , shown in the left column of 



of ~ 0, where the initial state is taken and where the 
absorption happens, has enormously small values being of 
the order of 10~^^. These amplitudes on the tail drop sig- 
nificantly with an increase of M. For the map ([3| the decay 
of eigenstates is more irregular since the band at x ~ 
crosses some secondary islands thus leading to appearance 
of a plateau in the decay with y. But in global we can still 
say that there is an exponential decay of eigenstates. This 
exponential localization of eigenstates reminds the Ander- 
son localization in disordered solid state systems (see e.g. 





10 20 . 30 40 50 



Fig. 15. (Color online) Modulus of the projection coefficients 
/jj of the initial density vector V^init, localized in one cell at 
xq = yo = 0.0625, with respect to the right eigenvectors ipj^ (of 
the UPFO projected case for M = 1600) versus level number 
j. These coefficients appear in the expansion V^init = fJ^j i^f 
(see text). The left and right panels represent data for the maps 
([1]) and ([3| respectively. The cases with — corre- 
spond to pairs of complex conjugated modes with /Xj+i = /i*. 



We can also consider the projection of our initial state 
taken in a cell ^0 on the eigenstates. Indeed, this initial 
state can be expressed as %l)\nit = l-^j "^f where fij are 
expansion amplitudes and iljf the right eigenvectors de- 
fined by Eq. ([2]). To determine the values of jij we need 
first to compute the left eigenvectors ip^ of the Ulam ma- 
trix Sp which are biorthogonal to the right eigenvectors 
ij)^ and provide the expansion amplitudes by the iden- 
tity: jij = <'0j"|'0init>/<'^j'|'^j^>- Note that this expres- 
sion does not depend on the choosen normalization of the 
eigenvectors and it requires only that < il)^\il)^> 7^ 0. 
However, for convenience, we have normalized both type of 
eigenvectors by the Li-norm such that ^ |?/^j^'^(x, y)\ = 
1. We have numerically determined the first 51 left eigen- 
vectors with the help of the Arnoldi method applied to the 
transpose of Sp and therefore obtained the corresponding 
expansion amplitudes. 



Fig. 14 , we see a clear evidence of exponential localization 
of eigenstates. In fact the average amplitude in a vicinity 



The dependence of 11 j on j is shown in Fig. 15 We see 
that there are enormously large fluctuations oi/ij which 
are in a range of 10 orders of magnitude. In particular 
the amplitudes corresponding to resonant modes are very 
small which is easy to understand if the resonant mode is 
localized far away from the initial state and does there- 
fore not contribute to the expansion. We think that these 
fluctuations are at the origin of the slow algebraic decay 
of Poincare recurrences P{t) (see below). 

In Fig. 16 we show the contribution of the largest Nm 
eigenmodes to the statistics of Poincare recurrences (for 
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Fig. 16. (Color online) Contributions of the largest eigen- 
modes of the UPFO projected case at M = 1600 to the statis- 
tics of Poincare recurrences for the maps ([T]) (left panel) and 
([3| (right panel). Here, we show the probability p{t) obtained 
from the expansion over eigenvectors given by the formula 
Pit) = Ef^o'ViAj with p, = M, E.,„^f(a;,y), Nm being 
the number of used modes and the eigenvalues being ordered 
as |Ao| > I All > IA2I > ••• (see text). The upper red curve 
is obtained from the direct iteration of the UPFO (see green 
curve in the right panels of figures [s] and and corresponds 
to the contribution of the full spectrum of all eigenvalues with 
Nm = Np. The middle blue curve corresponds to Nm = 51 with 
the same /ij values as those shown in figure [15] The main con- 
tributions to this curve arise from the diffusion modes (with 
real positive eigenvalues Xj > 0), the other resonant modes 
with complex or real negative eigenvalues give only a small 
contribution which does not modify the curve up to graphical 
precision. The bottom green curve corresponds to Nm = 1, 
i. e. the contribution fio Xq of the largest A eigenmode. In both 
panels the dashed line indicates for comparison a power law 
decay P(t) oc t"^ ^ 



M = 1600) given by the formula: p{t) = ^.^^^ ^ pj X] 
with Pj = iJLj y '4^f{x^ y) and the eigenvalues ordered 
as |Ao| > |Ai| > IA2I > 

For Nra = we have the statistics of Poincare re- 
currences obtained from the iteration of the UPFO and 
already shown in Figs. [S] and 12 For Nm = 51 we have 
evaluated the sum using the expansion coefficients shown 
in Fig. 15 Both curves coincide at t > 10^ for the map ([l]) 
or at t > 3 X 10^ for the map ^ showing that the largest 
eigenmodes determine the long time behavior. For large 
times {t > 10^ — 10^) only the first eigenmode contributes 
and the decay is purely exponential. It turns out that in 
the sum for Nm = 51 the terms arising from the resonant 
modes can be omitted without changing the curve up to 
graphical precision since these modes contribute only very 
weakly in the expansion. In general, the partial sum p{t) 
converges to the actual statistics of Poincare recurrences 
P{t) with increasing Nm and at given value of Nm one 
expects that p{t) and P{t) coincide for t 2 7^^ . 



The data of Figs. [Mj [15) [16] illustrate the nontrivial 
link between the localized eigenstates of the Ulam matrix 
and the decay of Poincare recurrences. The eigenmodes 
are exponentially localized and for many of them their 
projection on the initial state is very small but at some 
large times their contribution can become very important 
since the modes with large projections decay more rapidly. 



6 Discussion 

Our studies show that the generalized Ulam method re- 
produces well the decay of Poincare recurrences P{t) in 2D 
symplectic maps with divided phase space. At the same 
time the computation of P{t) is obtained in a more effi- 
cient way by the proposed SMCM allowing to reach time 
scales of the order of t = 10^^. We find that at these 
large times the Poincare exponent has values /3 = 1.58 for 
the Chirikov standard map at Kg and P = 1.70 for the 
separatrix map at Ac. The recurrences at large times are 
dominated by sticking of trajectories not only in a vicin- 
ity of the critical golden curve but also in a vicinity of 
secondary resonance structures. This confirms earlier nu- 
merical observations obtained on shorter time scales [27 . 

The sticking around various different resonant struc- 
tures on smaller and smaller scales of phase space leads 
to nontrivial oscillations of the Poincare exponent. The 
values of /3 found here are not so far from the average 
values found previously by averaging over maps at differ- 
ent parameters with ^ ^ 1.5 [18,19 , (3 ^ 1.57 [28,. In 
agreement with the data presented here and in [34], we 
find that the above value of f3 is close to the exponent of 
integrated density of states of the Ulam matrix which has 
P ^ 1.5. At the same time we see that at t = 10^^ the 
fluctuations in the Chirikov standard map at various Nf 
and various random realizations are significantly stronger 
as compared to the separatrix map. 

We attribute these fluctuations to a localization of 
eigenstates of the Ulam matrix which gives very nontriv- 
ial properties of eigenstates projection on an initial state. 
The properties of these eigenstates are still poorly under- 
stood. We think that the further developments of analyt- 
ical models of renormalization on Cayley type tree [22] 
28,32 ,33 and their applications to the puzzle of statistics 
of Poincare recurrences should develop a more detailed 
analysis of localization of eigenstates of the Ulam matrix. 
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